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Abstract Orbital motion of a body can be found from Newtonian equation of motion. 
However, it is useful to express the motion through time derivatives of Keplerian orbital 
elements, mainly if the motion is perturbed by small perturbing force. The first set of 
equations for the time derivatives of the orbital elements can be derived from the 
equation of motion using Lagrange brackets. The second one by using equation of 
motion and perturbation acceleration decomposed to radial, transversal and normal 
components. This paper shows that the second type of the perturbation equations can 
be derived from the first type using simple mathematical operations. 
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1 Introduction 

Perturbation equations of celestial mechanics belong to the standard part of the clas- 
sical astronomy. Perturbation acceleration represents disturbing force which disturbs 
Keplerian motion. Disturbing force causes time change of orbital elements describing 
actual orbit of a body in space. To describe the orbit of the body in space we can 
use different sets of orbital elements. We will use the following set of orbital elements: 
semi-major axis a, eccentricity e, inclination i, argument of perihelion ui, longitude 
of ascending node Q and angle a = nr, where n is mean motion and r is time of 
perihelion passage. Perturbation equations can be expressed using different methods. 
Using Lagrange method it is possible to derive perturbation equations which use scalar 
product of disturbing acceleration and partial derivative of position vector with respect 
to orbital elements. Another method is to derive perturbation equations directly from 
radial, transversal and normal components of disturbing force. We show that the ex- 
pression obtained by Lagrange method enables to derive the expression through radial, 
transversal and normal components of the disturbing force. Attempts to show the ex- 
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istence of this connection can be found in Brown (1896). Brown uses an alternate set 
of orbital elements and his equations contain several trivial errors. 



2 Expression obtained using Lagrange brackets 

Using Lagrange brackets, we can derive the following time derivatives of orbital ele- 
ments defined in previous section (see, e. g., Brouwer and Clemence 1961) 

da 2 dr 

777 = — & D ■ jt- , (1) 

at na do 

de 1 — e 2 dr \/l — e 2 dr . . 

-77 = o— ■ a o & D ■ TT~ , (2) 

at na z e do na z e ou> 

di cot i dr 1 dr . . 

dt = na 2yrr^2 ai3 ■ ~ na 2^Trr72 sin ~ az5 ' 9« ' ( ) 

do 2 dr \-e 2 dr 

-77 = 3-D ■ ~ o & D ■ 77- , (4) 

at na oa na z e oe 



duj _ Vl — e 2 <9r coti dr . . 

df2 1 dr 

= 2 /i 2 ~ a ^ ' 757 ' (6) 

at na^yl — e £ smj oi 

where ajj is a disturbing acceleration and r is a position vector of a particle with 
respect to the Sun. 



3 Expression through radial, transversal and normal component of 
disturbing acceleration 



We can express time derivatives of orbital elements though radial, transversal and 
normal components of disturbing acceleration in the following way 



dt 
di 
~dt 



da 
~dl 



de VT - ^ 



^V / ^ 



[a R e sin / + a T (1 + e cos /)] 



an sin f + ax [ cos / + 



e + cos / 
1 + e cos / 



= aN 



rcos(9 



na 

do _ 1 - e 2 
dt na 



2 V / T __ e 2 " 



aR 



cos/ 



1 + e cos / 



ax 



sin / 2 + e cos / 



e 1 + e cos / 



dn 
t— 
dt 



du> _ Vl — e 2 
dt nae 
dO 
HI 



—an cosf + ax sin/ 



r sin 



2 + e cos / 
1 + e cos / 



- a» 



r sin cot i 



a N 



na 2 Vl — e 2 sini 



(7) 
(8) 
(9) 
(10) 

(11) 
(12) 



where r = a(l — e 2 )/(l + ecos /) for an elliptical orbit and = f + u>. See e. g. Bate 
et al. (1971), Klacka (1992). 
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Fig. 1 Oscular orbital elements and radial, transversal and normal unit vectors. 



4 Relation 

When we want to evaluate time derivatives of orbital elements in Eqs. (l)-(6), we 
need to know the values of arj • (dr/dA), where A € {a, e, i, a, u>, O}. Moreover, we 
want to express time derivatives in Eqs. (l)-(6) through radial, transversal and normal 
components of perturbation acceleration arj = a R e R + a^e? + a N e N- To do this, we 
express also components of the vector u A = dr/dA through radial, transversal and 
normal components. We have 

dr dx . dy . dz . . 

UA = dA^dA 1+ dA i+ dA^^ W + MA « J + UAzk = 

= U AR e_R + u AT e T + u AN e N , (13) 

where i, j, k are unit vectors in directions of coordinate axes x, y, z of Cartesian 
coordinate system. Since vectors e^, ex, eN are orthonormal, we can write for radial, 
transversal and normal components of the vector u A 

uar = u Ax i ■ e R + u Ay i ■ e R + u Az k. ■ e R = u A ■ e R , (14) 
U AT = u Ax i ■ e T + u Ay i ■ e T + u Az \t ■ e T = u A ■ e T , (15) 
uan = u Ax i ■ e N + u Ay j ■ e N + u Az k ■ e N = u A ■ e N ■ (16) 

Finally, we get for arj • (dr/dA) 

dr 

a D ■ — = a R u AR + a T u AT + a N u AN ■ (17) 
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We need to find components of unit vectors e^, ex, ejy in Cartesian coordinate 
system as a function of orbital elements. We can write 

= (cosai,cosa2,cos«3) , (18) 

e T = (cos/?i,cos/?2,cos/?3) , (19) 

e N = (cos 71, cos 72, cos 73) , (20) 

where a\, a.2 and 0:3 are angles between vector eji and coordinate axes x, y and z, 
respectively. Similarly for vectors and e^r- To calculate one of the vectors e^, ey, 
e^v, we can use cross product of others two, since vectors are orthonormal. Components 
of the vectors can be found using a spherical law of cosines 

cos ZPOS = cos ZPOR cos ZROS + sin ZPOR sin ZROS cos S , (21) 

where each of the points P, R, S lie on one of the three different lines crossing in the 
point O (the point O is different from the points P, R, S) and S is the angle between 
planes determined by points P, O, R and R, O, S. 

We can find components of the unit vector e^j from spherical triangles AABD, 
ADBC and AEBD which can be constructed in Fig. 1. We will use notation — f +oj. 
From spherical triangle AABD we have 

ZBOA = n, ZDOB = 0, 6 ai = n - i, 

cosai = cos ZAOD= cos ZBOA cos ZDOB+ sin ZBOA sin ZDOB cos<5 Ql , 
cos o?i = cos Q cos — sin Q sin cos i, 

where S ai is angle between planes determined by points B, O, A and D, O, B. From 
spherical triangle ADBC 

ZDOB = 0, ZCOB = tt/2 - n, S a2 = i, 

cos a 2 = cos ZDOC= cos ZDOB cos ZCOB+ sin ZDOB sin ZCOB cos<5 Q2 , 
cos «2 = cos sin O + sin cos Q cos i, 

where 8 a2 is the angle between the planes determined by points D, O, B and C, O, B. 
From spherical triangle AEBD we obtain 

ZEOB = tt/2, ZDOB = 0, 6 a3 = tt/2 - i, 

cos a 3 = cos ZEOD= cos ZEOB cos ZDOB+ sin ZEOB sin ZDOB cos S a3 , 
cos«3 = sin sin i, 

where S a3 is the angle between the planes determined by points E, O, B and D, O, B. 

Components of unit vector ej\r can be found from spherical triangles AABF, ABCF 
and from the angle ZFOE depicted in Fig. 1. From spherical triangle AABF we have 

ZBOA = n, ZBOF = tt/2, <5 7l = tt/2 - i, 

cos 71 = cos Z AOE= cos ZBOA cos ZBOF+ sin ZBOA sin ZBOF cos <5 7l , 
cos 71 = sin J7sini, 

where <5 7l is the angle between the planes determined by points B, O, A and B, O, F. 
From spherical triangle ACFB we have 

ZBOF = tt/2, ZCOB = tt/2 - Q, <5 72 = tt/2 + i, 
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cos 72 = cosZFOC= cosZBOF cos ZCOB+ sin ZBOF sin Z COB cos5 72 , 
cos 72 = — cos Q sin i, 

where 5 72 is the angle between the planes determined by the points B, O, A and B, O, 
F. For ZFOE we have 

cos 73 = — cos ZFOE = cosi. 

Components of unit vector ex we can calculate using cross product ejv x eR- We 
can summarize the results as 

e R — (cos i? cos(/ + lo) — sin Q sin(/ + lo) cos i, 

sin J?cos(/ + ui) + cos J?sin(/ + lo) cosi, sin(/ + lo) sin i) , (22) 
ex = (— cos fi sin(/ + ui) — sin J? cos(/ + ui) cos i, 

— sin J?sin(/ + lo) + cos i7cos(/ + lo) cos i, cos(/ + oj) sin i) , (23) 
e N — (sin J7sin i, — cos sin i, cos i) , (24) 

where we have used = / + lo again. 

We begin with evaluating of da/dt. We need to calculate partial derivatives dx/da, 
dy/da and dz/da in Eq. (1). We can find components of the position vector from the 
relation r = r where 

_ a(l-e 2 ) 



1 + e cos / ' 

for the elliptical orbit. Finally, we get for x, y, z 



[cos ilcos(f + lo) — sin J?sin(/ + lo) cosi] o(l — e 2 ) 
1 + e cos / 



(25) 



(26) 



[ sin Q cos(/ + ui) + cos J? sin(/ + lo) cos i] o(l — e 2 ) 
1 + e cos / 



sin(/ + lo) sin i a(l- e 2 ) 



(28) 



1 + e cos / 

In these expressions only / is a function of a, thus we need to calculate df /da. For 
this purpose we use Kepler equation in the form 

nt + a = E-esmE , (29) 

where E is eccentric anomaly. We calculate partial derivative of Eq. (29) with respect 
to a and the result rewrite to the form 



dE 



da 1 — e cos E 

Now we use relation between the true anomaly and eccentric anomaly 



(30) 




/ = 2 arctan \ tan — . (31) 



Partial derivative of Eq. (31) with respect to a give 
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where we have used also the relation 



co S E=^L. (33) 
1 + e cos / 



We put Eq. (30) into Eq. (32) and the result we rewrite to the following form 

da r 2 

Now we can calculate partial derivatives dx/da, dy/da and dz/da 
dx 



(34) 



— Uax — 



| — cos fl [sin(/ + ui) + e sin w] 



da yi - ? 

— sin J? cos i [cos(/ + to) + ecoscj]} , (35) 

dy a r . . . . . , 

— = Uay — —=== \ — sin S2 |sm(/ + lu) + e sin u>\ 
da VI — e 2 

+ cos J? cos i [cos(/ + w) + ecoso;]} , (36) 

t^- = w<jz = a sini [cos(/ + w) + ecosw] . (37) 
da Vl — e 2 

When we put equations Eqs. (35)-(37) and Eq. (22) into Eq. (14), we get for the radial 
component of the vector u CT the following relation 

U<jR = "^j==f e Sin f ■ ( 38 ) 

Similarly, from Eq. (15) using Eqs. (35)-(37) and Eq. (23) we get for the transversal 
component 

u aT = jJL= (l + ecos/) , (39) 
VI - 

and, finally, for the normal component, we get, from Eq. (16) using Eqs. (35)-(37) and 
Eq. (24), 

UaN = ■ (40) 
Using Eqs. (38)-(40) in Eq. (17) we can now calculate a^> • (dr/da) 
dv a 

3l D ■ — = - [a R e sin / + a T (1 + e cos /)] . (41) 

Putting Eq. (41) into Eq. (1) we obtain for time derivative of the semimajor axis 

[a R esin/ + a T (l + ecos/)] . (42) 



This is the relation identical to Eq. (7). 

Now we want to calculate time derivative of eccentricity. We need to evaluate 
' (dr/duj) in Eq. (2). For partial derivatives dx/du>, dy/dui, dy/dui we obtain, from 

Eqs. (26)-(28), 

Q x \ — cos Q sin(/ + ui) — sin Q cos(/ + to) cos i\ a(l — e 2 ) 

■7T- = = — 7 , (43) 

duj 1 + e cos / 
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dy 



[ — sin i?sin(/ + to) + cos J?cos(/ + lo) cosi] a(l — e 2 ) 
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1 + e cos / 
cos(/ + u>) sini a(l — e 2 ) 



duj " 1 + e cos / 

This equations can be more simply written in our notation as (see also Eq. 23) 



dr 



a(l-e 2 ) 



duj 1 + e cos / 

Using Eq. (46) and Eqs. (14)- (16) we can immediately write 

UuR = , 
UojT = r , 
u^n = . 

By inserting Eqs. (47)-(49) into Eq. (17) we obtain 

dr 

a D ■ — = a T r . 



(44) 
(45) 

(46) 

(47) 

(48) 
(49) 

(50) 



Inserting Eq. (41) and Eq. (50) into Eq. (2) we obtain for the time derivative of the 
eccentricity 



de \/l — e 2 
dt na 



dR sin / + a T cos / + 



e + cos / 
1 + e cos / 



(51) 



This relation is identical with Eq. (8). 

Next we want to calculate di/dt. In this order we need to evaluate a^> ■ (dr/dfi) in 
Eq. (3). For partial derivatives of coordinates with respect to Q we get 



dx 



[ — sin J?cos(/ + to) — cos I2sin(/ + lo) cos i] a(l — e 2 ) 
1 + e cos / 

[ cos Q cos(/ + u>) — sin fl sin(/ + cu) cos i] a(l — e 2 ) 



1 + e cos / 



For radial, transversal and normal components of the vector uq we obtain 

U(2R = , 
uqt = r cos i , 
Uf2N — ~ r cos (/ + w ) sm * • 
By inserting Eqs. (55)-(57) into Eq. (17) we obtain 

&D ■ = r[ax cosi — cos(/ + u>) sini] . 

From Eq. (3), using Eq. (50) and Eq. (58), we get 

di rcos(f + Lu) 



(52) 

(53) 
(54) 

(55) 

(56) 
(57) 

(58) 
(59) 
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This relation is identical with Eq. (9). 

In order to calculate da/dt, we need to know sld ■ dr/da and ajj • dr/de in Eq. (4). 
We begin with a£> • dr/da. In Eqs. (26)-(28) we must take into account that also / is a 
function of a. Thus, we need to calculate df /da. We can use Eq. (31) to find similarly 
as in Eq. (32) 

!f=7CT (1+ecos/) S • (60) 

Now we use Kepler equation defined in Eq. (29). Partial derivative of Eq. (29) with 
respect to a gives the following relation 

dE _ 1 dn t _ 1 dn t 
da 1 — e cos E da 1 — e cos E da 



where n is the mean motion n = \J GM/a 3 (G is the gravitational constant and M is 
mass of central object). Using Eq. (34) (compare also Eqs. 32 and 30 with Eqs. 60 and 
61) we can write for df /da 

df dfdn^ . . 

da da da 

Using Eqs. (35)-(37) we can write for partial derivatives of coordinates x, y, z with 
respect to a 

dx [cosficos(/ + w) - sini?sin(/ + cj)cosi] (1 - (?) dx dn 

= Uax = — h -^--j-t , (63) 

da 1 + e cos / da da 

Qy [sin J?cos(/ + uj) + cos Q sin(/ + u>) cos i\ (1 — e 2 ) dy dn 

TT" = U ay = — h , (64) 

da 1 + e cos J da da 

dz sin(/ + lo) sini (1 — e 2 ) dz dn . . 

da 1 + e cos / da da 

These three equations can be more simple written in our notation as 

dr (i-e 2 ) , dn . ,aa\ 

tt" = "a = ^— 7 e R + - r tu tT ■ (66) 

da 1 + e cos / da 

By inserting Eq. (66) into Eqs. (14)-(16) we obtain 

(1 — e 2 ) dn ,„_. 

"afl = t~~ ^ + iVij -r-t , 67 

1 + e cos / da 

dn 

u a T = u aT -^t , (68) 

dn , . 

UaN = UaN -r^t . (69) 

When we put Eqs. (67)-(69) into Eq. (17) we get 

dr (1 — e 2 ) dr dn 

»d ■ 7T = «i? t 1 - + az, • tt- — t . 70 

da 1 + e cos / acr da 

We can now use Eq. (1) to obtain 

dr (1 — e 2 ) na da dn (1 — e 2 ) nadn. . 

azj • j- = a fl -i- + ___t = 0fl -i- ^7 + ^-^7*- (71) 

da 1 + e cos / 2 dt da 1 + e cos / 2 at 
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Moreover, we need to evaluate ajj ■ dr/de in order to find de/dt. In Eqs. (26)-(28) we 
need take into account that / is also a function of e. For partial derivative of Kepler 
equation defined by Eq. (29) with respect to e we obtain 



dE 
de 



sini? 



1 — e cos E 



(72) 



When we use Eq. (72) and also expressions a(l — ecosE) = r and <z\/I — e 2 smE = 
rsin/ in partial derivative of Eq. (31) with respect to e we finally obtain 



df_ _ sin/ 
de 



1 



(2 + ecos/) 



For partial derivatives of coordinates x, y, z we can, in our notation, write 



dr 



df 



de T de BT 



2ae 



o(l- 



(cos/-esin/|£) 



1 + ecos/ (1 + ecos/) 2 
This equation can be simplified using Eq. (73) as 



dr 



. 2 + e cos / 
a sin / ex — a cos / 



de 1 + ecos/ 

By inserting Eq. (75) into Eqs. (14)-(16) we obtain 

u e R = ~a cos / , 

. 2 + e cos / 

u eT = a sin / — - 

1 + e cos / 

u eN = . 

When we put Eqs. (76)-(78) into Eq. (17) we get 



dr 



—an a cos / + ax asm f 



2 + e cos / 



de 1 + e cos / 

Now we can use Eq. (71) and Eq. (79) in Eq. (4) to calculate da/dt we have 



da 



cos/ 



1 + e cos / 



ax 



sin / 2 + e cos / 
e 1 + e cos / 



- t- 



, dn 
~dt 



(73) 
(74) 

(75) 

(76) 

(77) 
(78) 

(79) 
(80) 



This relation is identical with Eq. (10). 

Now we want to calculate time derivative of argument of perihelion. In this order 
we need to evaluate &d ■ (dr/di) in Eq. (5). For partial derivatives dx/di, dy/di, dy/di 
we obtain from Eqs. (26)- (28) 



dx _ _ sin fl sin(/ + uj) sin i a(l — e 

"TTT — u ix — r. : 7 

dl 1 + e COS / 

dy _ _ — cos O sin(/ + lo) sin i a(l — e 

di ty 1 + e cos / 

dz _ _ sin(/ + u>) cos i a(l — e 2 ) 

di 1 " 1 + e cos / 



(81) 

(82) 
(83) 
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For radial, transversal and normal components of vector we obtain from Eqs. (14)- 
(16) 

u iR = , (84) 
u iT = , (85) 
u iN = rsin(/ + w) . (86) 
From Eq. (17), using Eqs. (84)-(86), we get 

Qy 

aD ' di = aN T Sin ^ + ^ ' ^ 
When we now put Eqs. (79) and (87) into Eq. (5), we obtain 

duj V 1 — e 2 / „ . r 2 + ecosf\ r sin( f + u>) cot i 

-„= -o« cos/ + a T sin/— L )-a N KJ , 88 

dt nae y 1 + ecos/y na 2 Vl - e 2 

This relation is identical with Eq. (11). 

Now we can finally put Eq. (87) into Eg. (6), in order to obtain dO/dt 

dO rsin(/+u) 

If ~ N 2 /i 2 ~ • (89 ' 1 

«i na^vl — e z sin! 

This relation is identical with Eq. (12). 



5 Conclusion 



We have just shown that it is possible to derive Eqs. (7)-(12) from Eqs. (l)-(6) using 
partial derivatives of position vector with respect to orbital elements. 
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